#library load
library(foreign)
library(base)
library(ggplot2)
library(weights)
library(reshape2)
library(survey)
library(mfx)
library(haven)
library(dplyr)
library(fmsb)
library(sf)
library(tidyverse)
library(sp)

#WD:
setwd() #set to your local working directory

#load shapefiles
#from https://data.humdata.org/dataset/mexican-administrative-level-0-country-1-estado-and-2-municipio-boundary-polygons
mx <- st_read("newshp/mex_admbnda_adm2_govmex.shp")

#load control municipalities
cont <- read.csv("control_municipalities.csv")

#load treatment municipalities
treat <- read.csv("treament_municipalities.csv")

#create control dummy:
cont$control <- 1
cont <-  cont %>% rename(ADM2_ES=municipio_residencia)

#create treatment dummy:
treat$treatment <- 1
treat <-  treat  %>% rename(ADM2_ES=municipio_residencia)


#merge control and treatment municipalities
dummies <- merge(cont,treat, by="ADM2_ES", all=TRUE) %>% mutate(control=if_else(is.na(control),0,1), treatment=if_else(is.na(treatment),0,1))

#correct case dummy Municipios
dummies$ADM2_ES <- str_to_lower(dummies$ADM2_ES)

#correct case dummy Municipios
mx$ADM2_ES <- str_to_lower(mx$ADM2_ES)


#correct municipio names
mx$ADM2_ES[mx$ADM2_ES=="chimalhuacrn"] <- "chimalhuacan"
mx$ADM2_ES[mx$ADM2_ES=="coacalco de berrioz"] <- "coacalco de berriozabal"
mx$ADM2_ES[mx$ADM2_ES=="coxcatlgn"] <- "coxcatlan"
mx$ADM2_ES[mx$ADM2_ES=="heroica ciudad de juchitln de zaragoza"] <-"heroica ciudad de juchitan de zaragoza"
mx$ADM2_ES[mx$ADM2_ES=="ju7rez"] <- "juarez"
mx$ADM2_ES[mx$ADM2_ES=="le0n"] <- "leon"  
mx$ADM2_ES[mx$ADM2_ES=="m5rida"] <- "merida"
mx$ADM2_ES[mx$ADM2_ES=="morole"] <- "moroleon"
mx$ADM2_ES[mx$ADM2_ES=="naucalpan de juirez"] <- "naucalpan de juarez"
mx$ADM2_ES[mx$ADM2_ES=="nicolts romero"] <- "nicolas romero"
mx$ADM2_ES[mx$ADM2_ES=="oaxaca de juorez"] <- "oaxaca de juarez"
mx$ADM2_ES[mx$ADM2_ES=="p2njamo"] <- "penjamo"
mx$ADM2_ES[mx$ADM2_ES=="rincjn de romos"] <- "rincon de romos"
mx$ADM2_ES[mx$ADM2_ES=="acapulco de ju rez"] <- "acapulco de juarez"
mx$ADM2_ES[mx$ADM2_ES=="almoloya de jurrez"] <- "almoloya de juarez"
mx$ADM2_ES[mx$ADM2_ES=="apatzingen"] <- "apatzingan"
mx$ADM2_ES[mx$ADM2_ES=="bah"] <- "bahaia de banderas"
mx$ADM2_ES[mx$ADM2_ES=="boca del rro"] <- "boca del rio"
mx$ADM2_ES[mx$ADM2_ES=="castaoos"] <- "castanos"
mx$ADM2_ES[mx$ADM2_ES=="r3o bravo"] <- "rio bravo"
mx$ADM2_ES[mx$ADM2_ES=="san andrss cholula"] <- "san andres cholula"
mx$ADM2_ES[mx$ADM2_ES=="san cristtbal de las casas"] <- "san cristobal de las casas"
mx$ADM2_ES[mx$ADM2_ES=="santa lucpa del camino"] <- "santa lucia del camino"
mx$ADM2_ES[mx$ADM2_ES=="texcaltitlcn"] <- "texcaltitlan"
mx$ADM2_ES[mx$ADM2_ES=="santo tomns hueyotlipan"] <- "santo tomas hueyotlipan"
mx$ADM2_ES[mx$ADM2_ES=="tihuatlbn"] <- "tihuatlan"
mx$ADM2_ES[mx$ADM2_ES=="santa maraa huatulco"] <- "santa maria huatulco"
mx$ADM2_ES[mx$ADM2_ES=="cuautitl n"] <- "cuautitlan"
mx$ADM2_ES[mx$ADM2_ES=="cuautitl n izcalli"] <-"cuautitlan izcalli"
mx$ADM2_ES[mx$ADM2_ES=="bahaia de banderas"] <-"bahia de banderas"
mx$ADM2_ES[mx$ADM2_ES=="cuauhtamoc"] <-"cuauhtemoc"
mx$ADM2_ES[mx$ADM2_ES=="san luis potosl"] <-"san luis potosi"
mx$ADM2_ES[mx$ADM2_ES=="c4rdoba"] <-"cordoba"  
mx$ADM2_ES[mx$ADM2_ES=="almoloya del rro"] <-"amoyala del rio"  
mx$ADM2_ES[mx$ADM2_ES=="atizapmn de zaragoza"] <-"atizapan de zaragoza"  
mx$ADM2_ES[mx$ADM2_ES=="atoyac de  lvarez"] <-"atoyac de alvarez"  
mx$ADM2_ES[mx$ADM2_ES=="c0rdenas" & mx$ADM1_ES=="SAN LUIS POTOSI"] <-"cardenas"  
mx$ADM2_ES[mx$ADM2_ES=="francisco i. madero"  & mx$ADM1_ES=="COAHUILA DE ZARAGOZA"] <-"francisco i madero"  
mx$ADM2_ES[mx$ADM2_ES=="cuauht moc"] <-"chuauhtemoc"  
mx$ADM2_ES[mx$ADM2_ES=="san juan bautista cuicatlhn"] <-"cuicatlan"  
mx$ADM2_ES[mx$ADM2_ES=="gustavo a madero"] <-"	gustavo a madero"  
mx$ADM2_ES[mx$ADM2_ES=="iguala de la independencia"] <-"iguala de la independencia"  
mx$ADM2_ES[mx$ADM2_ES=="l5zaro ccrdenas"] <-"lazaro cardenas"  
mx$ADM2_ES[mx$ADM2_ES=="minatitl n" & mx$ADM1_ES=="VERACRUZ DE IGNACIO DE LA LLAVE"] <- "minatitlan"  
mx$ADM2_ES[mx$ADM2_ES=="paratso"] <-"paraiso"  
mx$ADM2_ES[mx$ADM2_ES=="san luis r o colorado"] <-"san luis rio colorado"  
mx$ADM2_ES[mx$ADM2_ES=="san pedro garza garcea"] <-"san pedro garza garcia"  
mx$ADM2_ES[mx$ADM2_ES=="soledad de graciano stnchez"] <-"soledad de graciano sanchez"  
mx$ADM2_ES[mx$ADM2_ES=="tlajomulco de zueiga"] <-"tlajomulco de z??iga"  
mx$ADM2_ES[mx$ADM2_ES=="tonale"] <-"tonala"  
mx$ADM2_ES[mx$ADM2_ES=="torrean"] <-"torreon"  
mx$ADM2_ES[mx$ADM2_ES=="villa comaltitlln"] <-"villa comaltitlan"  
mx$ADM2_ES[mx$ADM2_ES=="benito jubrez" & mx$ADM1_ES=="QUINTANA ROO" ] <-"benito juarez" 
mx$ADM2_ES[mx$ADM2_ES=="benito juarez" & mx$ADM1_ES=="VERACRUZ DE IGNACIO DE LA LLAVE" ] <-"benito juarez1" 
mx$ADM2_ES[mx$ADM2_ES=="benito juarez" & mx$ADM1_ES=="ZACATECAS" ] <-"benito juarez1" 
mx$ADM2_ES[mx$ADM2_ES=="el arenal" & mx$ADM1_ES=="HIDALGO" ] <-"el arenal1" 
mx$ADM2_ES[mx$ADM2_ES=="iguala de la independencia"] <-"iguala de la independencia "  
mx$ADM2_ES[mx$ADM2_ES=="loreto" & mx$ADM1_ES=="BAJA CALIFORNIA SUR" ] <-"loreto1" 



#merge dummy to shp data
mx <- merge(mx,dummies,by="ADM2_ES", all=TRUE)


#add dummy for treatment
mx$control[mx$control==0] <- NA
mx$treatment[mx$treatment==0] <- NA
mx$type[mx$control==1] <- "Control"
mx$type[mx$treatment==1] <- "Treatment"

#municipality centroids
mx_pts <- st_centroid(mx)

#centroid coordinates
mx_pts <- cbind(mx_pts,st_coordinates(st_centroid(mx$geometry)))



#plot shp and centroids

#Figure 3:
pd <- position_dodge(0.1)
pdf("fg3.pdf", width = 10, height = 5.4)
map <- ggplot(fill="white") +
  geom_sf(data=mx, fill="white", color="grey70", size=.05) +
  geom_point(data=subset(mx_pts,control==1), aes(x=X,y=Y, color="Control"), size=2) +
  geom_point(data=subset(mx_pts,treatment==1), aes(x=X,y=Y, color="Treatment"), size=2)+
  scale_color_manual("Type", values=c(Treatment="darkred",Control="cornflowerblue"))+
  theme_void() + 
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()) +
  ylab("") +
  xlab("")
print(map)
dev.off()

